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We use spin-density-functional theory (SDFT) ab initio calculations to theoretically explore the 
possibility of achieving useful gate control over exchange coupling between cobalt clusters placed on 
a graphene sheet. By applying an electric field across supercells we demonstrate that the exchange 
interaction is strongly dependent on gate voltage, but find that it is also sensitive to the relative 
sublattice registration of the cobalt clusters. We use our results to discuss strategies for achieving 
strong and reproducible magneto-electric effects in graphene/transition-metal hybrid systems. 
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I. INTRODUCTION 

is an atomically thin two-dimensional 
gapless semi-conductor in which the carrier density can 
be varied over a broad range, from ~ — lO^'^ cm~^ to ^ 
-1-10^^ cm~^ by gating, and is a remarkably good conduc- 
tor at high carrier densities. Graphene/transition metal 
hybrid systems are attractive for spintronics because car- 
bon spin-orbit interactions are particularly weali^ in 
flat honeycomb-lattice arrays, because magnetic transi- 
tion element cluster^^^ form readily on graphene sur- 
faces, and because of potentially attractive propertie^^ 
of interfaces between graphene and magnetic transition 
metals. For example ultra-t hin tr ansition metal lay- 
ers on graphene are predictecP^'^ to have extremely 
large magnetic metal anisotropy energies. For these rea- 
sons there has recently been considerable interest^^^^ 
in the magnetic and electronic properties of transition 
metal adatoms and clusters placed on a two-dimensional 
graphene sheet. 

In this article we theoretically explore the possibility 
that the exchange coupling between separate magnetic 
metal clusters on graphene can be altered electrically by 
gating. Since arrays of magnetic clusters can be real- 
ized ongraphene by using a graphene/substrate moire 
patterrfl3 as a template, and the magnetic clusters hy- 
bridize relatively strongly with graphene's valence and 
conduction band orbitals, we anticipate gate-dependent 
exchange coupling between cluster s whic h should lead to 
gate-dependent magneto-resistanc^^^nil effects that are 
strong at room temperature. The goal of this work is 
to identify strategies for achieving strong, reproducible 
magneto-electric effects in graphene/transition-metal hy- 
brid systems. 

There is already a substantial theoretical 
literatur^20H28| Ruderman-Kittel-Kasuya-Yosida 
(RKKY) interactions between local moments coupled 
to graphene 7r-bands. It has been recognized,'22l for 
example, that when graphene is undoped the RKKY 
interaction is ferromagnetic (FM) for magnetic moments 
coupled to TT-electrons on the same graphene sublattice 
and antiferromagnetic (AFM) for moments coupled 



to TF-electrons on different sublattices. The RKKY 
interaction decays as at large distance r, because 
of the s uppressed density-of-states at the Dirac point of 
graphen6Pi23126128] At finite carrier density the RKKY 
coupling has spatial oscillations with period jhp on 
top of an envelope which decays as r~^. Most existing 
studies of the RKKY interactions in graphene have 
assumed magnetic moments due to point-like impurities 
that are associated with a particular honeycomb lattice 
site and have purely phenomenological interactions. 
These models are realized approximately in systems 
with magnetic moments due to hydrogenatiorP^ or 
carbon vacancie^^, although these defects significantly 
modify the carbon sp^ bonds and hence the structural 
and electronic properties of graphene. Moments due 
to adsorbed magnetic transition metal atoms do not 
distort the graphene bands as strongly but have small 
migration barrierJ^l! due to weak adsorption energies.l^ 
The transition metal clusters on graphene on which we 
focus are relatively immobile, however, and can be large 
enough to exceed the super paramagnetic limit. These 
larger magnetic objects therefore have more potential 
for spintronics applications. We attempt to realistically 
describe the magnitude of cobalt cluster moments, their 
magnetic anisotropy energies (MAE), the exchange 
coupling between the clusters and graphene, and finally 
the graphene-mediated magnetic exchange energies 
between separated clusters. 

We use first-principles supercell electronic structure 
calculations based on spin density functional theory 
(SDFT) to investigate not only the RKKY coupling be- 
tween magnetic cobalt clusters deposited on graphene, 
but also its dependence on external electric fields due to 
gating. We choose cobalt because its bulk lattice con- 
stant is very close to that of graphene, and because thin 
cobalt films down to two or three atomic layers have been 
found to have perpendicular magnetic anisotropj^^, which 
is preferable for spintronic applications. First, by calcu- 
lating the electronic structure of a two-atomic-layer thick 
cobalt film on graphene, we find that there is consider- 
able charge transfer from cobalt to graphene. Hybridiza- 
tion between the cobalt cluster and graphene leads to 
sublattice and spin dependent shifts in graphene 7r-band 
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energies from which we are able to extract the essen- 
tial kinetic-exchange parameters. Then we directly cal- 
culate the exchange interaction between two parallel two- 
atomic-layer-thick cobalt ribbons placed on graphene. 
For the geometries we have been able to consider, we 
find that the exchange interactions have a typical size 
~ 10""* eV per cobalt atom, comparable to the MAE of 
bulk cobalt (4 x 10~^ eVf^lt) ^nd thin films of cobalt on 
graphen^, but smaller than anisotropy e nergi es which 
can be achieved in asymmetrical clusters I^Mil g^jgQ 
find that exchange interaction tend to change sign when 
a graphene cluster changes its sublattice registration, and 
that the exchange interactions can be modified by gate 
voltages. 

In Section II we briefly describe the methods that we 
use for these computations. For the sake of dcfiniteness 
we have focused our attention on cobalt clusters that are 
two atomic layers thick and arranged in a ribbon geom- 
etry. In Section III we describe our results for the elec- 
tronic structure of a bulk two-layer thick film of cobalt 
on graphene. We find that there is considerable charge 
transfer from cobalt to graphene, and that hybridiza- 
tion between the magnetic cluster and graphene leads 
to sublattice and spin dependent shifts in graphene tt- 
band energies. In Section IV we summarize our results 
for the dependence of total energy on the relative spin 
orientations of separated clusters. We are able to under- 
stand our main findings using an approximate treatment 
which treats the cobalt-graphene interaction pcrturba- 
tively. Finally in Section V we present our results for the 
gate-voltage dependence of these exchange interactions. 
We find that gate fields can produce sizable changes in 
exchange interactions, in some cases changing their signs 
and substantially reducing their sublattice registration 
dependence. In Section VI we summarize our findings 
and discuss some possible directions for future research. 



culations. Denser fc-point meshes (up to 79 x 79 x 1) were 
used to check accuracy and to perform MAE calculations. 

To study the indirect exchange coupling between re- 
mote cobalt clusters on graphene, we constructed a su- 
percell with two parallel cobalt ribbons two-atomic-layers 
thick and three atoms wide, oriented along the zigzag di- 
rection of graphene (Fig. [2]) . The supercell used in this 
case is 25 X 1 with the same 20 A vacuum layer in z direc- 
tion. These ribbon calculations used a 1 x 49 x 1 fc-point 
mesh. The lattice parameters of the cobalt ribbons were 
taken from the infinite 2D slab calculations mentioned 
above without further relaxation. (We checked the in- 
fiuence of relaxation for several cases and did not find 
qualitative modification relative to the results reported 
on below.) The exchange coupling between the cobalt 
ribbons was estimated by calculating the total energy dif- 
ference between spin-parallel (EM) and spin-antiparallel 
(AFM) configurations: 

Ai? = il^FM — -EafM- (1) 

With this convention a positive Ai? corresponds to anti- 
ferromagnetic exchange between the ribbons. 

An external electric field across the supercells was re- 
alized by adding a saw-tooth like external potential to 
the total energy functionaF^. We have applied electric 
fields of different size in the same supercell as in Fig. [2] 
In this case the external field can produce only charge 
transfer between the two cobalt ribbons and graphene. 
A more realistic representation of gating action on the 
graphene/transition metal hybrid system can be achieved 
by adding a bilaycr Cu slab to the supercell as in (Fig. |9| . 
The copper acts as a a charge reservoir and also screens 
the part of graphene directly below the cobalt ribbons 
from external fields. A more detailed discussion of some 
issues involved in using VASP to simulate gates is pro- 
vided in Appendix B. 



II. METHODS 

The DET calculations reported on in this work were 
performed using the projector-augmented-wave (PAW:P 
method as implemented in the Vienna ah initio simula- 
tion package (VASPj^^^^. The Perdew-Burke-Emzerhof 
generalized gradient approximation (PBE-GGA)"^^ was 
used for the exchange-correlation energy functional. To 
calculate the electronic band structure of an infinite 
graphene sheet fully covered by a two-atomic-layer-thick 
cobalt film [Fig. [l] (a)], we used a a 20 A thick vacuum 
region between neighboring supercells in the z (perpen- 
dicular to the graphene plane) direction. We fixed the 
lattice constant at the experimental value for graphene 
(2.46 A) since the (0001) surface of bulk hep cobalt has a 
small lattice mismatch (< 2%). All atoms in the supercell 
were allowed to relax until the Hellmann-Feynman force 
on each atom was smaller than 0.001 eV/A. A plane- wave 
energy cutoff of 400 eV and a 33 x 33 x 1 fc-point mesh 
were used for structure relaxation and total energy cal- 



III. KINETIC EXCHANGE COUPLING 
BETWEEN COBALT OVERLAYERS AND 
GRAPHENE tt-BANDS 

A. Ah Initio Spin-density-functional Theory 

As illustrated in Fig. [T] (a) , we have calculated the to- 
tal energies of bilayer cobalt films adsorbed on graphene 
with different registries and have found that the most 
stable geometry is that with the C atoms in one sub- 
lattice of graphene located directly below bottom-layer 
cobalt atoms, i.e. at atop sites, and the C atoms in the 
other sublattice below the top-layer cobalt atoms, i.e. at 
hep sites. The optimal separation between the cobalt 
overlayer and graphene is about 2.21 A. After adsorb- 
tion on graphene, the magnetic moments on the cobalt 
atoms in the first layer (adjacent to graphene) decrease 
from 1.710 per cobalt atom, which is close to the bulk 
value, to 1.560 /ub per cobalt atom. Meanwhile, the C 
atoms in sublattice A (adjacent to cobalt atoms) obtain 
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a per-atom magnetic moment of 0.043 ^b, antiparallel 
to the magnetization of the cobalt overlayer, whereas the 
C atoms in sublattice B acquire a moment of 0.041 /iB 
per atom and parallel to the cobalt moments. There- 
fore the overall magnetization direction of graphene is 
opposite to that of the cobalt film. We have also calcu- 
lated the magnetocrystalline part of the MAE by eval- 
uating the total energy difference, including spin-orbit 
coupling, between configurations with all moments along 
the z direction (out-of-plane) and along the x direction 
(in-plane). The system is found to have perpendicular 
magnetic anisotropj^^, with a MAE of ~0.09 meV per 
cobalt atom, which is larger than that of bulk hep cobalt 
(~0.04 meV), but still the same order of magnitude. 

The spin-resolved Kohn-Sham band structure of the 
Co-graphene hybrid system is shown in Fig. [l] (b). The 
graphene bands are spin-split and the Dirac points at 
the K point are gapped because of the relatively strong 
interaction with the cobalt overlayer, in agreement with 
previous result^^^''^. It is nevertheless clear from the 
position of the Fermi level that graphene is n-doped, 
i.e. electrons are transferred from cobalt to graphene^. 
The graphene layer majority-spin Dirac point is easily 
identified in the two-dimensional bands, but its minor- 
ity spin-counterpart is so strongly hybridized with cobalt 
d-orbitals that it is less easily identified. The K point 
is at a higher energy for graphene majority spin bands 
than for minority spin bands, indicating an overall anti- 
ferromagnetic coupling between the cobalt overlayer and 
graphene. This conclusion is also in agreement with the 
antiparallel orientations of the graphene and cobalt mag- 
netizations mentioned above. 
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FIG. 1: (color online), (a) Top and side views of the su- 
percell (with a 3 x 3 repetition in the xy plane for illustra- 
tion purpose). The larger balls represent cobalt atoms and 
the smaller balls C atoms, (b) Two-dimensional Kohn-Sham 
quasiparticle band structure of the Co-graphene hybrid sys- 
tem neglecting spin-orbit interactions. The blue lines illus- 
trate the majority spin bands and the red lines the minority 
spin bands. The blue and red dots indicate the strength of 
carbon orbital character in the majority and minority spin 
states, (c) Model graphene projected band structure calcu- 
lated using Eq. [2] The model parameter values (Eq.JsJl are 
obtained by fitting to the DFT results listed in Table lir 



B. Kinetic Exchange Model 

Our electronic structure calculations can be qualita- 
tively described using a simple model for graphene cou- 
pled to a cobalt overlayer in which hybridization and 
charge transfer effects shift the energies of both majority 
and minority spins on both graphene sublattices: 

H = hvpk - T + ^i- ho.zTz - h^fiSz - h^^:,SzTz. (2) 

In Eq. [2]the first term on the right hand side is the usual 
Dirac Hamiltonian for hopping on a honeycomb lattice 
with velocity vp ^ 10® m/s and wave vectors measured 
relative to the Brillouin-zone corners, Sz = ±1/2 labels 
spin, and Tz — ±1 distinguishes A (under the atop site) 
and B (under the hep site) sublattices. The parameters 
of this model can be identified by fitting to the energies 
of the bands that have the largest 7r-band character at 
the Brillouin-zone corner (K) points, which are summa- 
rized in Table |lj is diagonal when fc = and its four 
eigenvalues 
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The four Kohn-Sham bands with the strongest carbon 
Pz character at the K point of Brillouin zone are bands 
2, 3, 4, and 5 in Table |Tj By fitting their energies to 
the SDFT band energies we can obtain the values of the 
parameters: 

fi ^ -0.622 eV (5) 
ho,z = 0.195 eV 
hz,o = -0.214 eV 
hz.z = -0.766 eV. 

The model band structure calculated with these these 
parameters is plotted in Fig. [l] (c) . 
Several comments are in order: 

(i) The chemical potential fi specifies the energy shift av- 
eraged over spin and sublattice, which is negative because 
electrons are transferred to graphene, in agreement with 
our previous discussion. 

(ii) The value of /ig^ is positive because the A sublattice 
is more strongly influenced by the cobalt overlayer than 
the B sublattice, which is expected since the A sublattice 
is directly below the cobalt atoms at the interface. 

(iii) The value of hz^ measures the kinetic exchange 
coupling between cobalt and graphene spins averaged 



4 



over sublattices. Its negative sign means the sublattice- 
averaged magnetic coupling is AFM, also in agreement 
with our observations in the previous subsection. 

(iv) The spin- and sublattice-dependent term hz^z reflects 
the property that the majority spin is higher in energy 
on the A sublattice whereas the minority spin is higher 
in energy on the B sublattice. In other words, the Co- 
graphene exchange coupling is AFM on the A sublattice 
but FM on the B sublattice. To understand this property 
we refer to Table|l]in which bands 1, 6, 7, 8 are identified 
as cobalt d bands that hybridize with the graphene tt 
bands. From the carbon pz and cobalt d characters that 
these bands carry, it can be seen that spin-splitting on 
the A sublattice is because of hybridization mainly with 
the d"^^ orbitals of cobalt (bands 1 and 6), whose 
minority spin states are higher in energy than majority 
spin states and above the Fermi level. The higher energy 
of the carbon majority spin states on the A sublattice can 
therefore be understood as the result of level repulsion 
from cobalt d"^^ orbitals with the same spin. The 
same argument also applies for the B sublattice, whose 
Pz orbitals mainly hybridize with the d^^, d^^ , d^^, and 

2 2 

d^ ~y orbitals of cobalt because of symmetry. However, 
both of the two cobalt d bands (band 7 and 8) with these 
characterare below the Fermi energy and the 7r-bands 
at the K point, with the minority spin band higher in 
energy. Therefore level repulsion in this case results in 
the higher energy of the minority spin states, i.e. in 
ferromagnetic coupling. 

(v) hz^z is much larger than hz^ because the kinetic ex- 
change interaction between the cobalt overlayer and the 
graphene is strongly dependent on sublattice. We will 
see later that this property will translate to a strong de- 
pendence of the graphene-mediated exchange interaction 
between two cobalt clusters on their relative registries 
with respect to the sublattices of a continuous graphene 
sheet. 



TABLE I: Orbital character of the bands in Fig. [T] (b) at 
the K point of 2D Brillouin zone. Only those having strong 
carbon characters are listed. A and B correspond to the 
two sublattices of graphene, as shown in Fig. [l] (a) . 



Band No. Energy (eV) C Pz character cobalt d character 


1 


1.006 


Ai: 0.069 


3^2 - i: 0.671 


2 


-0.151 


Bi: 0.319 


xz, yz i: 0.098 






x^ - J/^a;y i: 0.099 


3 


-0.328 


At: 0.297 


6z — r t- 0.3DO 


4 


-0.703 


Bt: 0.439 


xz, yz t: 0.059 






x^ - y'^,xy t: 0.016 


5 


-1.307 


a;: 0.341 


3z^ - i: 0.019 


6 


-1.754 


At: 0.191 


3z^ - t: 0.303 


7 


-1.965 


Bi: 0.185 


xz, yz i: 0.19 






x^ - y'^,xy i: 0.054 


8 


-3.048 


Bt: 0.055 


xz,yzt: 0.206 






x^ - y'^,xy t: 0.166 



FIG. 2: (color onhne). Top and side views of the supercell (re- 
peated by 4 times in the y direction for illustration purposes) 
used to calculate the magnetic couphng between two parallel 
cobalt ribbons (larger blue balls) placed on a graphene sheet 
(smaller yellow balls). 



A. Electronic Structure 



IV. MAGNETIC COUPLING BETWEEN 
COBALT CLUSTERS ON NEUTRAL GRAPHENE 



In this section we will investigate the magnetic cou- 
pling between cobalt clusters on neutral graphene sheets 
which are mediated mainly by their mutual influence on 
the graphene 7r-bands. First we employ SDFT to study 
a relatively small system with parallel quasi-lD cobalt 
ribbons placed on graphene (Fig. [2) and separated by 
~ 1 nm. Then we will calculate the RKKY coupling 
in graphene perturbatively using the the model devel- 
oped above to compare with the SDFT calculation re- 
sults. This comparison informs perturbative estimates 
of coupling which cannot be directly addressed using ab 
initio tools. 



In Fig. [3] (a) we show the electrostatic potential (ionic 
potential plus Hartree potential from electrons) profile 
within the graphene sheet for the system in Fig. [2] In 
equilibrium, the chemical potential will shift relative to 
the bands by the opposite amount. Therefore Fig. [3](a), 
with a sign change and up to a constant, can be viewed as 
a plot of Ferm energy relative to the Dirac point. One can 
see that there is a large positive shift of chemical poten- 
tial in the region directly below the two cobalt ribbons, 
meaning the graphene is strongly n-doped at these posi- 
tions. The TT-band electron barrier height between cobalt- 
covered and bare graphene regions is therefore about 0.5 
eV, close to the 0.622 eV separation between the chemical 
potential and the Dirac point found earlier for the infinite 
2D Co/graphene hybrid system. The barrier is smaller in 
the present case because separations between neighboring 



5 



cobalt ribbons are not large enough for the pristine neu- 
tral graphene value. This barrier can potentially decrease 
magnetic coupling between remote graphene clusters by 
localizing electronic states more strongly in the vicinity 
of one particular cluster. 

In Figs. [3] (b-d) we plot partial density-of-states 
(PDOS) functions projected to the Pz orbitals of car- 
bon atoms at different points in the structure. At all 
three sites the PDOS Dirac-point minima are shifted to 
lower energy, indicating n-type doping over the entire 
graphene sheet. The magnitude of the Dirac-point shift 
decreases as one goes further away from the cobalt rib- 
bons, as expected. One feature worth mentioning in the 
PDOS plots is the appearance of resonant features that 
are absent in pristine graphene. These features can be 
identified as confinement effects in the zigzag-ribbon-like 
uncovered graphene regions between the cobalt ribbons. 
We see later that although these modifications to the lin- 
ear DOS of graphene do not greatly influence the form of 
the TT-band mediated magnetic coupling, they do play a 
role in the dependence of the charge transfer to graphene 
on gate field. 
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FIG. 3: (color online), (a) Electrostatic potential variation 
in adsorbed cobalt ribbons, (b-d) Density of states projected 
to the Pz orbitals of three carbon atoms whose positions are 
indicated by the black arrows. Black lines-graphene with ad- 
sorbed cobalt ribbons, red lines-bare graphene. The negative 
PDOS axis plots minority band values while the positive axis 
plots majority band values. 



B. Exchange Coupling 

We next study the exchange coupling between the two 
cobalt ribbons in Fig. [2] In Fig. |4] we plot the spin 
density vs. position within the graphene sheet for the 
case of two ferromagnetically aligned cobalt ribbons. In 
the region below the cobalt ribbons, the spin polariza- 
tions are opposite for the two sublattices of graphene, as 



in the case of complete two-layer cobalt coverage. This 
property is maintained in the uncovered portion of the 
graphene sheet. Opposite spin polarizations on the two 
sublattices suggests that the graphene-mediated inter- 
action will be strongly sublattice dependent as in the 
RKKY case. Th is be havior is common in systems with 
bipartite lattices.l^^ISll 
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FIG. 4: (color online). Color scale plot of spin polarization as 
a function of position within the graphene plane, in the region 
between two cobalt ribbons with parallel spin orientations. 
The vertical axis in this figure is on position along the ribbon 
direction which has atomic scale periodicity. The positive and 
negative spin densities (in arbitrary units) are concentrated on 
carbon atoms on opposite sublattices. The black dots indicate 
the positions of C atoms. 

In Fig. [5] we plot SDFT results for magnetic coupling 
between cobalt ribbons for different edge-to-edge separa- 
tions between the ribbons and different registries with 
respect to the sublattices of the continuous graphene 
sheets. We first note that although both cobalt ribbons 
have the same atop-hcp registry with graphene, the first 
layer cobalt atom is sometimes atop an A site carbon 
atom and sometimes atop a B site carbon atom. The con- 
figurations of atop(A)-hcp(B) and atop(B)-hcp(A) are 
degenerate for an individual cobalt ribbon, but magnetic 
coupling energies can change if one ribbon changes reg- 
istry and the other does not. The strong oscillation be- 
tween FM and AFM coupling in Fig. [s] is due to pre- 
cisely this effect. From now on we refer to the geometry 
in which the two cobalt ribbons have the same registry 
or different registries respectively as geometry AA, and 
geometry AB. 

From Fig. [5] we see that the strength of the magnetic 
coupling is about 1.3 meV per supercell for the AA con- 
figuration for separations between 8 A and 17 A. This 
exchange coupling is about 0.13 meV when normalized 
per cobalt atom, which is much larger than the 0.04 meV 
MAE of bulk hep cobalt and somewhat larger than the 
MAE of a 2 layer cobalt film on graphene (0.09 meV). 
(We have also calculated the MAE of a single cobalt rib- 
bon on graphene as in the present setup and the value 
is 0.08 meV per cobalt atom, with the easy axis along 
the ribbon direction.) The similar strength of the MAE 
and the exchange coupling means that inter-ribbon inter- 
actions can have a substantial influence on the magnetic 
conflguration of cluster arrays. RKKY-like oscillations in 
the coupling are expected to have period ^ vr/Zci?, with 
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kp the Fermi wave vector. In the present system the 
Fermi energy Ep is about 0.4 eV on average in the part 
of graphene between the two cobalt ribbons, correspond- 
ing to a period of ^ 5 nm. Therefore it is not surprising 
that we do not see RKKY-hkc oscihations in these calcu- 
lations. The small coupling at distances below 5 A may 
be due to competition between direct exchange coupling 
and graphene-mediated coupling between the two cobalt 
ribbons. It is not clear why there is strong variation 
in the exchange coupling strength for the AB configura- 
tion. One guess is that it is due to structural details at 
the boundaries of the zigzag-ribbon-like graphene region 
between the two cobalt ribbons. 
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FIG. 5: (color online). Magnetic coupling (per supercell, 
which has 10 Co atoms) between two cobalt ribbons as a 
function of ribbon separation. The interaction strength is 
the total energy difference between parallel and antiparallel 
spin- alignment configurations. Black squares (red dots) cor- 
respond to configurations in which the cobalt atoms in the 
bottom layers of the two cobalt ribbons are directly above 
the same (different) sublattice(s) of graphene. 

It is important for potential applications to understand 
how these exchange couplings will change with the size 
of the cobalt clusters. Due to computational power lim- 
itations we consider only two cases. First we increase 
the width of the two cobalt ribbons from 3 to 4 atoms, 
so that there are 14 cobalt atoms in a supercell. In the 
second case we add one more layer of cobalt atoms to 
the 4-atom-wide ribbons in case 1, so that the number of 
total cobalt atoms increases to 18. The per-cobalt mag- 
netic coupling is 0.10 and 0.099 meV for the two cases. 
In both cases the per-atom coupling strength is smaller 



than the 0.13 meV value obtained at the original cluster 
size. Therefore one can expect the total exchange cou- 
pling to increase sub linearly with cluster size. There are 
several reasons why this finding is expected. First, as we 
mentioned previously, there is a large chemical potential 
barrier at the cluster edge, which will weaken the influ- 
ence of cobalt atoms deeper inside the clusters. Second, 
when the cluster size is comparable to or larger than the 
oscillation period of the RKKY interaction, contributions 
from different parts of the cluster interfere destructively, 
as we see in the next subsection. Finally, since the largest 
contribution to the kinetic exchange interaction between 
cobalt clusters and graphene is from the cobalt atoms 
closest to graphene, adding more layers of cobalt to the 
clusters is expected to be less effective in increasing the 
magnetic coupling. 



C. Qualitative Theory of Exchange Coupling 

In this subsection we will use conventional perturba- 
tion theory and the model defined by Eq. [2] to calcu- 
late the RKKY coupling between magnetic clusters on 
graphene, and compare the result with our SDFT re- 
sults. Similar calculations for the RKKY interaction in 
graphene has been performed previously,!^^}!!!] mainly 
for the case of point-like magnetic impurities. Here we 
will explicitly include the size and shape of magnetic clus- 
ters. When combined with the essential kinetic exchange 
parameters obtained from first principles, the formalism 
developed in this subsection can be a useful tool for ex- 
trapolations to system sizes beyond the range which cov- 
ered by SDFT calculations. 

For a graphene sheet that is partially covered by two 
distinct magnetic clusters 1 and 2, Eq. [2] becomes 

H = H0 + H1+H2 (6) 
= hvFk-T + Di{r)Yi+D2ir)Y2, 

where 151(2) (r) = 1 at positions covered by cluster 1 
(2) and zero otherwise, and Vi(2) = /^i(2) ^ ''■o.z'''z,i(2) ~ 
hz,oSz^i{2)-hz,zSz,i{2)Tz,ii2)- Therefore the RKKY inter- 
action is evaluated by calculating the contribution to the 
total energy at second order in the perturbation H1+H2: 



I 



where g — 2\s the valley degeneracy, s — ±1 is the band 
index, and fgk is the Fermi distribution function [1 -I- 
exp((i?sfe — ^) /kBT)]~^ . In keeping with the continuum 
model we are using to describe the graphene 7r-bands, we 



neglect inter-valley transitions which add an anisotropic 
and rapid modul ation to the spatial dependence of the 
RKKY interactioii2I128l. 
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The eigenfunctions of Hq are: 

1 /g-iffc 



in which Dg ijo) is the Fourier transform of £>i(2)(''') 
Therefore Eq. [T] becomes 



ikr 



(8) 



where 0^ = arctan(fcj,/fcj;). ^^1(2) can be written as a 
Fourier integrah 



Vl(2)Vl(2) 



(9) 



{2n) 



2 ifsk fs'k-\-qj 



Esk — Es'k+q 



(10) 



By substituting Eq. |8] and the spin-dependent terms in 



V 



1(2) 



into |Fj,fc^^(7^q,iVi +-Dq,2V2)F,fcr, and keeping 

'i and Dq^2^2, we 



only the cross terms between Dq,!^ 
obtain 



l^^i'fc + qPq.lVl + i?q,2V2)F,fc|2 = (11) 
iD*q,lDq.2 + C.C.) • {^hl oil + SS' COs{0k - Ok+q)] 
+ 7;hl.z[^ ~ Ss' C0a{9k - Ok+q)]T^^iT^^2}S:,.lSz.2, 



AE 



in which the first term in the curly brackets is 
sublattice- independent and the second term is sublattice- 
dependent. Here 1(2) are ±1 depending on which 
graphene sublattices the clusters are directly above. For 
conciseness we set Sz,iSz,2 1/4 from now on. Note 
that the cross terms between hzfi and hz.z vanish be- 
cause unperturbed graphene has spatial inversion sym- 
metry, and Tz changes sign under spatial inversion. Us- 
ing the values of /iz,o and hz_z obtained previously, the where 



factor multiplying the sublattice-dependent term is ~ 10 
times larger than that the factor which multiplies the 
sublattice-independent term. Therefore the RKKY in- 
teraction between cobalt clusters should be strongly de- 
pendent on their registration with respect to the sublat- 
tices of graphene, agreeing with our observation from the 
SOFT resuhs. 

The integration over k and the summation over bands 
in Eq. [TO] can be performed explicitly at T = K. 
(We summarize calculation details in Appendix A.) The 
RKKY energy, written as an integral over q, is 



(2) 

B.KKY 



gh 



2 



\Qhvp 



(2^)2 



(i?;ii?q^2 + c.c.)n,,o(di2) 



iDq^Dq,2 + C.C.)Ilz,z{q)Tz.lTz,2, 



n.,o('7) 



kp 



kp 
2^ 



A 



kp 
vr 



27r 



2kp 

q 

2kp 

q 



2kf 



2kf 



e{q~2kp) + ^e{2kp-q), 



arcsin Q{q — 2kp) — -Q{2kp — q), 



(13) 
(14) 



Q{x) is the Heaviside step function, and A is the Dirac 
model's ultraviolet cutoff. Note that both 0(9) and its 
first derivative are continuous at q = 2kp. In contrast, 
Hz^ziq) has a discontinuous first derivative at q = 2kp, 
similar to the behavior of 2-dimensional electron gas. 
Therefore one can expect that the contribution to the 
RKKY interaction from the sublattice-independent part 



will have a faster decay with distance than that from the 
sublattice-dependent part.l^^ 



Graphene's RKKY interaction can be obtained by set- 
ting Di{r) = 6{r) and 1)2 (r) = S{r - R). The kpR > 1 
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limit is 



IRKKY 



(R) 



2 

z.z 



i?3 



when graphene is undoped, and 



J. 



RKKY 



(R) 



ghl^kp sm{2kpR) 



i?2 



,(15) 



(16) 



when graphene is doped. When carriers are present the 
dominant contribution is the sublattice-dependent part, 
which is oscillatory in space and decays as R^^. When 
graphene is undoped, the oscillatory term vanishes be- 
cause of the kp prefactor, and the leading order terms 
monotonically decay as R~^. 



Next we use Eq. 12 to calculate the RKKY-like inter- 



action between two cobalt ribbons on graphene with the 



same geometry as in our SDFT calculations. The distri- 
bution functions for this case are 

A(r) = e(^x + ^+0e(^-:r-0 (17) 

D2{r) = Q(^-x + w+^^e(^x-^ 

where d is the distance between the inner edges of the 
two ribbons, and w is the width of the two ribbons. Their 
Fourier transforms are 



Dg.l = 
Dq,2 = 

Therefore, 



1 

1 



2TrS{qy). (18) 
2T:S{qy). (19) 



^q,l^q,2 + C.C. = ^ {2 COs[qx {d 
^x 



— cos{qxd) — cos[qx{d + 2w)]} ■ 2nLS{qy). 



(20) 



In deriving the above equation we have used the relation 



S'iqy) = S{0)S{qy) 



L 
2tt 



(21) 



where L is the length of the system in y direction. We 



can then carry out the integration in Eq. 12 numerically. 
Below we will compare the results from this model cal- 
culations to the SDFT results. Note that the interaction 
energy in SDFT is the difference between spin-parallel 
and spin-antiparallel configurations of the two cobalt rib- 
bons. Therefore the model results below are all double 
Erkky inEq.[l2] 

Fig. [6] shows the magnetic coupling from our model 
for the AA geometry, which correspond to Tz^iTz,2 = 1 
in Eq. [T2j One can see that the order of magnitude 
agrees very well with the SDFT resuhs in Fig. [5] and 
the trend with changing distance is also well reproduced. 
We have chosen Ep to be 0.4 eV, which is the average of 
the graphene chemical potential under the cobalt ribbons 
(^0.6 eV) and that in the center between the two cobalt 
ribbons (~0.2 eV). The agreement would be improved if 
we sued the fact that the doping level of the graphene re- 
gion between the two cobalt ribbons increases as the two 
ribbons approach to each other. In Fig.[6](b) we assumed 
simple linear dependence of Ep with d and the agreement 
with Fig. [5] is remarkably improved. We note here that 
there is some arbitrariness in determining the width of 
the cobalt ribbons w since there is no sharp boundary of 
the portion of the graphene region which interacts with 
the cobalt ribbon. Here we chose to be 4 unit cells 
of graphene to account for the residue influence at the 
edges of the cobalt ribbons, although in a pure geomet- 



rical sense the cobalt ribbon amounts to 3 unit cells of 
graphene. w may be treated as a fitting parameter in 
applications of our approximate theory. The magnetic 
coupling for the AB geometry {tz^iTz^2 = —1 in Eg. [l2|, 
which we did not show in Fig. |6j can be obtained simply 
by subtracting the sublattice-dependent part from the 
sublattice-independent part. As we mentioned before, 
the anomalous oscillation of magnetic coupling for the 
AB geometry in Fig. [5] probably has a structural origin 
that is not captured by this simple model. 

Knowing that our model can capture the essential 
physics of the graphene-mediated magnetic coupling be- 
tween cobalt clusters relatively well, we can now explore 
the large separation limit which cannot be easily ad- 
dressed by first-principles methods. First in Fig.[7](a) we 
plot the magnetic coupling for the AA geometry vs. rib- 
bon separation for several carrier densities. One can now 
see the spatial oscillation between AFM and FM inter- 
actions which appears only beyond the separation range 
covered in Fig. [5j From the figure we see that not only 
the periodicity, but also the amplitude of the oscillation, 
depends on the doping level. This behavior is consistent 
with the asymptotic RKKY interaction Eq. [16] 

In Fig. [7] (b) we plot magnetic coupling divided by rib- 
bon width w, which is proportional to the magnetic cou- 
pling per cobalt atom. It is interesting to see that when 
w is very large (24 graphene unit cells in the zigzag direc- 
tion, equivalent to about 50 A), the magnetic coupling 
is strongly suppressed. This behavior can be understood 
by considering in terms of destructive superposition be- 
tween different parts of the ribbon, when the scale of 
the clusters is close to the oscillation period. In addi- 
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FIG. 6: (color online) (a) Model results for RKKY-like cou- 
pling between cobalt ribbons. _Bf=0.4 eV, hvF=5.9Q eV-A, 
L=2.46 A, m;=8.51 A. (b) Same as (a) but with Ep increasing 
linearly as d decreases. 



tion, since the period of the RKKY oscillation increases 
with decreasing kp, the coupling for the same large clus- 
ters will be less suppressed as the graphene is less doped, 
which we have also verified. Fig. [t] (b) also confirms our 
discussion on the effectiveness of increasing the magnetic 
coupling by preparing larger clusters. Therefore a gen- 
eral criterion for real applications is that the linear size 
of the clusters should be around or below which is 
half of the RKKY period. 



V. GATE CONTROL OF EXCHANGE 
COUPLING 

Since the RKKY coupling in graphene has a strong 
dependence on the Fermi energy (Eq. [16] and Fig. [?]), 
which in turn can be altered by electric gates, we expect 
that the magnetic coupling between cobalt clusters can 
be conveniently tuned by gating. In this section we will 
study the change of the magnetic coupling between cobalt 
ribbons on graphene with external electric fields. We 
have relegated some general remarks on how to simulate 
electric gates in supercell calculations to Appendix B. 



A. Freestanding Co-graphene in an Electric Field 

By directly applying a electric field along the z direc- 
tion in the supercell of Fig. [2j we can change the Fermi 
energy in the graphene by transferring electrons from the 
cobalt ribbons to graphene and vice versa. In Fig. [8] we 
show the charge transfer within the supercell after ap- 
plying a 0.2 V/A electric field along the —z direction. It 



FIG. 7: (color online) (a) RKKY coupling between cobalt 
ribbons at large separations, for several carrier densities, (b) 
RKKY coupling divided by ribbon width w for several widths. 
w is expressed in terms of the number of graphene unit cells 
along the zigzag direction across the cobalt ribbon, w is fixed 
at 4 in (a) and Ep is fixed at 0.4 eV in (b). 



can be seen that electrons are transferred from graphene 
to Co, and that an out-of-planc polarization is induced 
in the graphene sheet itself. The amount of charge trans- 
ferred from the graphene plane decreases as one moves 
away from the cobalt ribbons, in agreement with the elec- 
trostatic potential profile shown in Fig.[3](a). In this way 
one decreases the graphene carrier density not only in the 
bare regions of graphene, but also in the regions covered 
by the cobalt ribbons. 

One question which may be raised at this point is 
whether or not the exchange coupling between cobalt 
and graphene will be infiuenced by the electric field. 
To this end we have calculated the spin polarization in 
a graphene sheet fully covered by a 2-layer cobalt film 
[Fig[l](a)] under electric fields up to 0.8 V/A and did not 
find a significant change. Therefore the field dependence 
of the exchange coupling between graphene and cobalt 
is not an issue in the range of electric fields considered 
here. 

Next we study the field dependence of the magnetic 
coupling between the two cobalt ribbons at specific sep- 
arations between them. In Fig. [9] (a) we plot magnetic 
coupling vs. electric field for two cobalt ribbons sepa- 
rated by ^15 A, and different registries with the graphene 
sublattices. One can see that both the sign and magni- 
tude of the magnetic coupling can be tuned by electric 
fields. It is also interesting to notice that for both the 
AA and AB configurations the coupling has a similar 
sublinear dependence on electric field. Using the sim- 
ple model explained in the previous section, we found 
that the coupling changes almost linearly with Ep from 
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Electric field (V/A) Energy (eV) 



FIG. 8: (color online). Charge density difference (in an x — z 
plane) between a system subjected to a 0.2 V/A electric field 
along the —z direction, and a system with no electric field. 
Positive and negative values (in arbitrary units) correspond 
to accumulation and depletion of charge, respectively. The 
black dots (triangles) indicate the positions of C (Co) atoms 
in the plane. 



Ep — 0.2 eV to 0.4 eV, vifhich is roughly the range of 
Ep shift produced by the electric fields in our DFT cal- 
culations [Fig. |9] (b)]. Therefore the nonlinearity should 
come from the field dependence of the Fermi energy of 
graphene. In equilibrium the external potential difference 
between cobalt and graphene (eEd where d is the spatial 
separation) should be balanced by the electric potential 
due to charge redistribution and the Fermi energy shift 
of graphene (i.e., the quantum capacitance of graphene). 
This screening physics can be described crudely using a 
simple parallel plate capacitor model: 



ed-dE^ 



edcEp ■ dEp 
C 



+ dEf 



(22) 



lin- 



where c is the proportionality constant for the 
ear dependence of graphene DOS on Ep, and c = 
27r(tof )2 =0-018 eV~^A~^ in pure graphene, C/d is the 
geometric capacitance of the graphene/cobalt bilayer, 
and dE and dEp are electric field and Fermi energy dif- 
ferentials. The solution of this differential equation is 



Ep 



V2e2cd2C ■E + C^ + 2ecd ■ const - C 
ecd ' 



(23) 



which explains the slower-than-linear dependence Ep 
on E. Of course this argument relies on the assumption 
that the density of states of graphene around Ep is linear 
in energy. By looking at Fig. [o] (b) one can see that this 
assumption is actually reasonable, although the effective 
value of c may be different from that in pure graphene 
value due to the confinement-induced resonances. 



Finally in Fig. 10 (a) we plot magnetic coupling 



the separation between the two cobalt ribbons for several 
electric field strengths. The corresponding result from 
the model in Sec. IV C is plotted in Fig. [lO](b). Reason- 
able agreement for the E ~ —0.4 V/A case is obtained 
by taking Ep = 0.36 eV, which means this extremely 
large electric field is only able to shift Ep by 0.04 eV on 
average. The small number is partly due to the incom- 
plete coverage of the cobalt ribbons on graphene, which 



FIG. 9: (color online), (a) Dependence of magnetic coupling 
between two cobalt ribbons on external electric field at two 
different separations. Blue squares (red dots) correspond to 
the configuration that the two cobalt wires sit above the same 
(different) graphene sublattice(s), with a separation of 15.0 A 
(14.3 A). A negative value of field strength means that the 
field is along the ~z direction, (b) Density of states (spin- 
up plus spin-down) projected to the pz orbital of a C atom 
in the center of the supercell for several different external 
electric field strengths and the AA configuration in (a). The 
inset blows up the details around Ep. 



decreases the effective capacitance, but mostly due to the 
small vertical separation between the two systems, which 
makes graphene's quantum capacitance effect dominant. 
It is clear that an external electric field does not ade- 
quately model the infiuence of a remote gate. In the 
next subsection we will use an alternative supercell to 
better simulate a realistic gating geometry, and find that 
this tactic brings additional benefits. 



> 

E 




Separation (A) 



FIG. 10: (color online), (a) Magnetic coupling between two 
cobalt ribbons in the AA configuration vs. separation, under 
different electric fields, (b) Results obtained using the model 
in Sec. HVCl 
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B. Co/graphene with a Cu Slab Mimicking a Gate 
Electrode 



Fig. [TT] shows an alternative supercell which simulates 
electric gating more realistically. A two atomic layer 
thick slab of Cu is inserted in the supercell, at a distance 
of about 4 A from the graphene sheet. Because of its high 
density of states, the Cu slab will act as an electron reser- 
voir, just like a real gate electrode. We apply the electric 
field on the cobalt side of the graphene sheet and place 
the Cu slab on the other side of the sheet. Electrons are 
then transferred to or from the bare regions of graphene 
from the Cu slab, depending on the sign of the electric 
field. The part of graphene sheet that is directly below 
the cobalt ribbons is shielded from the the electric field 
by cobalt-layer screening. Consequently, complications 
due to field-dependent graphene cobalt coupling are mit- 
igated. Our calculations were motivated by the expecta- 
tion that adding carriers to the uncovered portion of the 
graphene sheet would reduce the potential barrier at the 
cobalt ribbon edges and in this way enhance magnetic 
coupling. 




FIG. 11: (color online). Top and side views of the supercell 
with a bilayer Cu slab (grey balls) mimicking a backgate. 
The supercell is repeated four times in the y direction for 
visualization purposes. 



In Fig. 



12 we show PDOS for different C atoms in 
the graphene sheet when no external magnetic field is 
applied. By comparing with Fig. [3] one can see that 
the PDOS is changed mainly by a shift of ~0.1 eV to- 
wards higher energies, which means that graphene is less 
n-doped. This result may seem counterintuitive since 
graphene is also n-doped on Cu, and Cu has an even 
smaller work function than that of Co. However, the di- 
rection of charge transfer when separation exceeds the 
range of direct chemical interaction is determined by rel- 
ative work functions. Because Cu has a larger work func- 
tion than graphene, it p-dopes graphene when chemically 
isolatecP2l. The p-doping by Cu enables us to explore a 
doping range of graphene that cannot be easily reached 
by directly applying an electric field to the freestanding 
Co-graphene system as in the previous subsection. 

Fig. [l3l (a) shows the charge transfer after applying 
a 0.2 V/A electric field along —z direction. One can 
see that electrons are indeed transferred from the Cu 
slab to the graphene and cobalt system. The part of 
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FIG. 12: (color online), (a-c) Density of states projected to 
the pz orbitals of three carbon atoms, for the supercell with 
a Cu slab. Black lines-graphene with the cobalt ribbons on 
top and the Cu slab below, red lines-bare graphene. 



graphene directly below the cobalt ribbons has almost 
no charge transfer, whereas the bare regions of graphene 
are electron-doped. The overall effect is essentially the 
same as would be produced by gating action from a pla- 
nar electrode separated vertically by a distance smaller 
than the graphene ribbon width. From the electrostatic 
potential plot in Fig. 13 (b), the potential barriers in 
graphene due to the cobalt ribbons are indeed reduced 
after applying the field ('-^O.OB eV by aligning the poten- 
tial at the cental region). The change is small because 
much of the external field is screened by the cobalt rib- 
bons and the Cu slab. This is a limit set by our supercell 
size, and is therefore an artifact of our calculation pro- 
cedures, but cannot be easily circumvented. Screening 
of the gate field due by metal clusters on graphene will 
however, be important experimentally when the distance 
to the gate is larger than the cluster separation. 

Because of the different charge transfer behavior in the 
present supercell compared to that without the Cu slab. 
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the field dependence of the magnetic coupling [Fig 
(a)] is changed. Without applying the electric field, the 
magnetic coupling is reduced because of the lower carrier 
density in the graphene between the two cobalt ribbons, 
as we have discussed previously. However, when a 0.2 
V/A field is applied along the —z direction, the coupling- 
separation curve is changed by reduced barrier heights. 
Namely, when the barrier height is lower, the increase 
in the average doping in between the two cobalt ribbons 
when they get closer will be less dramatic. Since the 
coupling is roughly proportional io kp^ the shape of the 
coupling-distance curve should be more tilted to the left. 
The scenario is consistent with the model explained in 
Sec. UVCl 

On the other hand, when a 0.2 V/A field is applied 
along the z direction, the coupling-distance curve is rela- 
tively smooth below 13 A, a behavior which we are able 
to reproduce using our model. A large shift of the curve 
appears at around 14 A. A tentative explanation is the 
following: When the distance between the two cobalt rib- 
bons is large, the central graphene region between them 
is nearly neutral. (In Fig. 14 (b) we show the PDOS 



of a carbon atom at the central region between the two 
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FIG. 13: (color online), (a) Charge density difference be- 
tween the systems subjected to a 0.2 V/A electric field along 
the —z direction, and no electric field. Positive and negative 
values (in arbitrary unit) mean accumulation and depletion 
of charge, respectively. Black dots, triangles, and squares in- 
dicate the positions of C, Co, and Cu atoms in the plane, 
respectively, (b) Relative electrostatic potential as defined in 
Fig. [3] (a), for systems subjected to a -0.2 V/A electric field 
(red lines), and zero electric field (black lines), respectively. 



cobalt ribbons, and it is seen that the DOS is almost hn- 
ear with energy.) Therefore Eq. 23 also apphes, according 



to which the change of Ep with field will be more pro- 
nounced when Ep is small. This effect, together with 
the fact that graphene will be more exposed to the exter- 
nal field as the two cobalt ribbons move away from each 
other, will likely lead to a sudden change of magnetic 
coupling at a certain separation. 




Separation (A) 



Energy (eV) 



FIG. 14: (color online), (a) Magnetic coupling vs. separation 
between cobalt ribbons in the AA configuration, for several 
electric fields. Fields are in units of V/A. (b) Density of states 
projected to the orbital of a C atom in the center of the 
supercell, with and without external electric fields, for the AA 
configuration and the separation of 17.1 A. The inset blows 
up details around Ep. 



VI. DISCUSSION AND CONCLUSIONS 

In this study we have demonstrated that cobalt mag- 
netic clusters on graphene can have relatively strong 
gate-voltage-dependent exchange interactions, but that 
these interactions are sensitive to the relative sublattice 
registration of cobalt clusters with respect to a mono- 
lithic graphene honeycomb. Although we have focused on 
cobalt clusters, the combined SDFT and phenomenologi- 
cal modeling approach used here can be straightforwardly 
applied to other systems, e.g. Ni clusters on graphene. 
We have carried out some similar calculations for Ni clus- 
ters, and find they have weaker exchange coupling with 
graphene than Co clusters. Thus cobalt has the distinct 
advantages of having both large exchange coupling and 
a good lattice match with graphene. 

In Fig. |3]and Fig. [12] we have seen that resonances in 
the density-of-states appear due to the quantum well and 
edge states of the zigzag-ribbon-like uncovered graphene 
segments in our supercell calculations. Although these 
density-of-states resonances do not have overwhelming 
importance for exchange interactions in the parameter 
range we were able to explore in this work, the phenom- 
ena may be interesting in their own right. For exam- 
ple, it is known that ideal graphene zigzag ribbons have 
spin-polarized edge states,H3nS31 but that graphene rib- 
bons with impurity-free edges are very difficult, if not 
entirely impossible, to fabricate experimentally-^SZI The 
study of spin-polarized graphene edge states, resulting 
from parallel magnetic ribbons deposited on graphene, 
may be an alternative route to realizing the potentially 
interesting edge physics of graphene nanoribbons. 

Our study addressed only the case of atop-hcp registry 
of the cobalt clusters with respect to graphene. This 
is the structure assumed by large 2D cobalt films on 
graphene. Given the strong sublattice registration de- 
pendence of this interface structure, we anticipate sim- 
ilar sensitivity to other structural modifications. We 
conclude that for any nanoparticle assembly method, 
precise control of the interface structure, at least in 
the first atomic layers, will be a crucial issue if repro- 
duceable exchange interactions are desired. In particu- 
lar, cobalt nanoparticles prepared using wet chemistry 
method^S^HUl are not likely to have consistent interface 
structures, and are therefore likely to have highly vari- 
able interactions. Her e we note that some authors have 
concluded theoreticalljl^'^ that the atop- fee interface be- 
tween graphene and Co(OOOl) is energetically slightly 
preferred to atop-hcp. The difference relative to our cal- 
culations could be due to a cobalt film thickness depen- 
dence of the preferred registry, or even due to differences 
in the exchange-correlation potentials used in the DFT 
calcualtions. Nevertheless, we found that the graphene- 
Co exchange coupling for the atop-fcc configuration does 
not differ qualitatively from the atop-hcp configuration, 
which is expected since the dominant contribution to the 
exchange coupling between the Co clusters and graphene 
is from the atop surface Co atoms. The structural ar- 
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rangement of the first row of magnetic atoms is however 
crucial. 

From our calculations, we can identify several key 
parameters that will influence the experimental real- 
ization of interesting magneto-resistance and magneto- 
electric devices in graphene/magnetic- metal hybrid sys- 
tems. Ideally we would like to be able to substan- 
tially alter the magnetic configuration of a cluster ar- 
ray by changing a gate voltage. For this to happen, 
the inter-cluster exchange coupling should be strongly 
gate- voltage-dependent and the same order of magnitude 
as the MAE. For clusters of fixed shape, we can expect 
that the per-atom MAE (~ 10~* eV) should be roughly 
cluster-size independent. The per-atom exchange cou- 
pling depends on cluster size, inter-cluster distance, and 
gate voltages. We can conclude that per-atom exchange 
coupling will be comparable to the MAE only for rela- 
tively small cluster sizes, and for relatively small inter- 
cluster distances. A reasonable bound for the inter- 
cluster distance is the period of the RKKY oscillation 
Tr/fep, which is on the order of a few nm for graphene 
with a large carrier density. The cluster size also must be 
smaller than this number to avoid destructive superposi- 
tion of coupling from different parts of a cluster. There- 
fore, the system size considered in our SDFT calculations 
is actually close to the ideal scale for strong effects. This 
length scale is obviously difficult to achieve, and will lead 
to magnetic and magneto-electric hysteresis only below 
^ lOO-R'. As we mentioned in the introduction, graphene 
moire patterns on metal substrates provide one attrac- 
tive strategy to achieve patterning on this length scale. 
These systems would have the disadvantage, however, 
that there would be no control over the relative sub lat- 
tice registration between different cobalt clusters. An- 
other strategy is to grow large domain graphene sheets on 
cobalt thin films and then etch away the metal connect- 
ing different regions. In this case it should be possible to 
maintain control over relative sub lattice orientation, but 
reaching the required length scales would be challenging. 

It is interesting to compare the related case of interac- 
tions between magneti c clu sters mediated by topological 
insulator surface states.'^^'^In both cases the 2D metallic 
states are described by a Dirac model. The main differ- 
ences in the topological insulator case are that graphene's 
sublattice degree of freedom is absent and that spin- 
orbit interactions are strong. Both differences point to 
potential advantages of the topological insulator struc- 
tures. The strong spin-orbit interactions at the TI sur- 
face will lead to strong magnetic anisotropies both in the 
energies of individual magnetic cluster^^, and in their 
interactional^, which will assist hysteresis at smaller clus- 
ter sizes. Most importantly, the absence of a sublattice 
degree of freedom should make it easier to control the 
magnetic interactions between clusters. 

In summary, we have described a survey of graphene- 
mediated exchange coupling between cobalt magnetic 
clusters, and of its tunability via electric gates. Our anal- 
ysis is based on ab initio SDFT calculations interpreted 



using approximate models. By fitting SDFT calculations 
of the electronic structure of a 2D thin film of cobalt de- 
posited on a single layer graphene sheet to a phenomeno- 
logical kinetic exchange model, we have identified the rel- 
evant kinetic exchange coupling parameters. From these 
parameters we were able to establish that the exchange 
coupling between cobalt clusters is strongly sublattice 
registration dependent. We then directly calculated the 
magnetic coupling between two infinite long cobalt rib- 
bons on graphene using SDFT, and found that their cou- 
pling is of the same order as the magnetic anisotropy 
energy of the cobalt ribbons. As expected, the coupling 
is found to change dramatically as one changes the rela- 
tive registries of the two cobalt ribbons with the graphene 
sublattices. We also identified the large potential barrier 
at the edge of the cobalt ribbons, which may influence the 
magnetic coupling in a variety of ways. To explore the be- 
haviors of the magnetic coupling in a much larger param- 
eter range, we constructed a phenomenological theory of 
the magnetic coupling using the simple Dirac Hamilto- 
nian of graphene and the kinetic exchange parameter we 
had obtained from the 2D calculations. The RKKY cou- 
pling given by this theory agrees well with the DFT re- 
sults for the same system. We found that the magnitude 
of the coupling depends on the Fermi energy of graphene, 
and that the coupling per cobalt atom will actually be 
very small when the cluster size is very large. By ap- 
plying an electric field inside the supercell in our SDFT 
calculations, we found that the electric field can lead to 
a considerable change in both magnitude and sign of the 
magnetic coupling between cobalt ribbons. The coupling 
changes faster with field when graphene is less doped, 
which was explained as a capacitance effect. We were 
also able to use the phenomenological theory to capture 
these behaviors. To better simulate the realistic gating 
configuration, we put a Cu slab in the supercell mim- 
icking a backgate, which also suppresses the potential 
barrier at the edge of the cobalt ribbons. We found that 
the change of coupling with field becomes more sensitive 
to the separation between the two ribbons, which is a 
consequence of the reduced potential barriers. 
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Appendix A: RKKY coupling of graphene from the 
continuous Dirac model 



We start from calculating the integral in Eq. [TO) taking the sublattice-dependent term as example: 



{Dl,Dg^2 + C.C.) 



d^fc 



(27r)2^ 

We first consider the situation of T = K and Ep at the Dirac point. Define 



,r r 1 - SS' C0s{9k - 0k + q) 

2 [Jsk — Js'k+q) 77771 _,|7. , -j Tz^iTz^2- 



s\k\ - s'\k + q\ 



'_(rv _ rO . 1 - Ss'cOs{9k - Ok+q) 

(2^^2U«fc Js'k+q) s\k\-s'\k + q\ ' 
where |(1 — s). To evaluate this integral we will have to calculate fl" ^{q) at finite freqency ui: 



1 - Ss'c0s(6'fe - Ok+q) 



_ fO X 

■'''''+''^ s\k\- s'\k + q\+uj + iS 



E 



dkde 



rUk- 



k+q cos t 
\k+q\ 



(27r)2 uj + a{k + \k + q\)+iS' 
where S is a small real number, a = ±1, and then take the limit of a; — >■ Ol^The result is 

nL(9)-|-A, 

where A is a cutoff. Similarly, for the sublattice- independent part, we got 



n°,o(9) 



which agrees with previous result^^lIini6T]_ 

Next we consider the doped case. Still take the sublattice-dependent part as example, and let 



fs'k+q) 



1 - SS' C0s(6>fc - 9k+q) 

s\k\ - s'|fc + q\ 



= 2E 



d fc ~ 1 - ss' cos(6'fc - 9k+q) 

Jsk 



s\k\ - s'\k + q\ 



(Al) 



(A2) 



(A3) 



(A4) 



(A5) 



(A6) 



dkd9 ^ kcosO 
(27r)2-^'= 9 + 2A:cos6l 

where fsk = fsk — fsk' fk = fii^ik) + fi{Eik + 2/i). The integral can be done straightforwardly. The result is 

(A7) 

For the sublattice-indepedent part, after similar calculations, we got 
kp kp 



An^ M ^ ^ _ A arcsin —<d{q - 2kp) - ^e(2fcF - q). 
n 2tt q 4 



An,^o(9) = 



TT 27r 



'1 - 



2kp 

q 



q . 2kp 

— — arcsin 

2kp q 



e{q~2kp) + ^e{2kp ~ q). 



(A8) 



Next we study the behavior of the graphene RKKY interaction between two point defects, namely, the RKKY 
range function. The distribution function is now Di{r) = S{r) and D2{r) = S{r — R), and 



Dl^Dg,2 + Dl^Dq^i = 2cos{q ■ R). 



(A9) 
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First we consider the sublattice-independent part. When graphene is undoped, i.e., kp 
Therefore we have 



0, n,,o('?) 



(Eq.lASi 



IGhvp 
16hvF 



d^q q cos(q • R) 



(27r)2 
dq q^ 



4 

MqR): 



(AlO) 



where Jq is the 0th order Bessel function. To evaluate this integral we refer to the formulaP^ 
The result is 



gh 



z,0 



1 



(All) 



(A12) 



Therefore at zero doping the sublattice-indepedent part corresponds to an antiferromagnetic interaction, and goes 
like at large R. 

The situation is a little complicated when graphene is doped. Since 0(9) is not singular at q = 2kp, the asymptotic 
behavior of Jz.o{R) at large R should be largely determined by the value of Ilzfi{q) at small q. However, Ilzfi{q) = 
when q < 2kp (Eq. A5 and A8). Therefore we can argue that Jz.o{R) in the doped case is a superposition of two 



terms with similar magnitude. However, one term (corresponding to the kp=0 contribution) decays monotonically as 
R~^ without oscillation, while the other term will be oscillating with the periodicity related to kp since there will be 
singularity at g = 2A:i;- in the higher order derivatives of Ilz.o{q). Therefore the long range behavior of Jzfi{R) should 
still be approximately proportional to and modulated with some oscillation. 
Finally we turn to the sublattice-dependent part Jz,z{R) 



Jz,z{R) 



gh 



2 

z,z 



dq 
27r 



2qIlz,z{q)Jo{qR)Tz,lTz^2- 



Note that 



n.,.(g) 



TT 



q . '^kp 

— arcsm fc)(g — 2kp) 

1-K q 



\Q{2kp~q) 



(A13) 



(A14) 



2 2kp , ^, 

1 arcsm kdiq — 2kp) 

vr q ' 



where we have dropped the constant terms since their Fourier transform will just be delta functions centered at i? = 0. 
The result of the integral is expressed in terms of the Meijer G-function: 



gK 



2 

z.z 



1 



-G: 



3,0 
2,4 



1,1 

3 3 1 



{kpRf 



(A15) 



16/it)i.7rii?3 ^'*V0'2'2'2 

The asymptotic behavior of Meijer G- functions at large argument can be found, e.g., in Ref. i62i We finally obtain 
the asymptotic form of Eq. |A15| at kpR^ 1 

ghl 



'^z,z 



1 



16hvp n^R^ 



cos{2kpR) + kpRsm{2kpR) 



(A16) 



The asymptotic expression of the Meijer-G function turns out to work very well (Fig. 15 ). So the sublattice-dependent 



contribution to the RKKY interaction has an oscillating form w ith th e period n/kp, and the leading order term decays 
as R~^, similar to the behavior of two-dimensional electron gaJ23EZl. 

We can finally write down the expression for the RKKY range function in graphene at kpR 3> 1 by keeping only 
the leading order term: 



J, 



RKKY 



(R) = 



ghl^kp am{2kpR) 
IGn^fwp R^ 



Jrkky{R) 



gh 



z,0 



1 



gh 



Tz.iTz^2 (doped), 
1 



T^,iT2,2(undoped). 



(A17) 
(A18) 



r 



Appendix B: Simulating gates in supercell 
calculations 



cell calculations. For this purpose it is natural to assume 



In this appendix we briefly discuss some of the chal- 
lenges in realistically simulating gates using VASP super- 
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- Meijer G-function 

- Asymptotic formula 




FIG. 15: Meijer G-function in Eg. [AIS] and its asymptotic formula in Eg. |A16| 



a slab geometryl^SEll. An external potential in the super- 
cell can be modeled straightforwardly by adding its inter- 
action energy with electrons and ions to the Kohn-Sham 
energy functional. However, because the potential corre- 
sponding to a homogeneous electric field is unbounded in 
space, to recover the periodic boundary condition of the 
supercells one needs to compensate the potential differ- 
ence between neighboring supercells. The usual scheme 
to do this is to add a fictitious dipole layer in the vac- 
uum, at the boundary of the supercelP^. The size of the 
dipole should be determined self-consistently in the min- 
imization process of the Kohn-Sham functional, so that 
the dipole layer will compensate the jump of the total po- 
tential rather than that of the external potential alone. 
The dipole layer must be homogeneous laterally, so that 
it will not induce artificial fields applied to the system 
of interest inside the supercell. As a result, shifting the 
system in the supercell as a whole towards or away from 
the dipole layer should in principle have no impact on 
the properties of the system itself. In other words, the 
external field in the supercell is like that from two gates 
at plus and minus infinity, respectively. 

This feature is not desirable when one would like to 
simulate a circumstance in which a laterally inhomoge- 
nous system, like our graphene sheets partially covered 
by cobalt ribbons, that is close to a gate. The surface of 



a real gate is an equipotential surface, so that charge will 
redistribute on it when the gate is close to a system that is 
laterally inhomogeneous. One strategy to simulate such 
an equipotential boundary condition is to place a a real 
metal slab inside the supercell. However, attention must 
be paid to another difference between supercell DFT cal- 
culations and real gates, i.e. that all subsystems share 
the same chemical potential in the VASP case. This is a 
result of energy minimization in solving the Kohn-Sham 
equation by taking the whole supercell as one system. 
Consequently, spatially separate parts in the supercell act 
as if they were all electrically shorted. We have utilized 
this property in Sec. |VB| In the slab geometry we con- 
sidered here, anything between two metal slabs (provided 
that they are thick enough) will be screened from exter- 
nal fields. Therefore, the best choice to simulate a real 
gate close to a system is to shift the system close to one 
boundary of the supercell, and put the metal slab at the 
opposite boundary from the system. We have tried this 
geometry using the supercells considered in this paper 
and found it indeed works well. The geometry, however, 
will not do better than the supercells used in the main 
text, in terms of the simulating cases with large shifts in 
graphene Fermi. This is because the charge redistribu- 
tion on the metal slab will actually decrease the field felt 
by the regions of graphene not covered by cobalt ribbons. 
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